home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / src / rand.cc < prev    next >
C/C++ Source or Header  |  1997-07-19  |  9KB  |  413 lines
  1. /*
  2.  
  3. Copyright (C) 1996 John W. Eaton
  4.  
  5. This file is part of Octave.
  6.  
  7. Octave is free software; you can redistribute it and/or modify it
  8. under the terms of the GNU General Public License as published by the
  9. Free Software Foundation; either version 2, or (at your option) any
  10. later version.
  11.  
  12. Octave is distributed in the hope that it will be useful, but WITHOUT
  13. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  15. for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with Octave; see the file COPYING.  If not, write to the Free
  19. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  20.  
  21. */
  22.  
  23. /* Modified by Klaus Gebhardt, 1997 */
  24.  
  25. #ifdef HAVE_CONFIG_H
  26. #include <config.h>
  27. #endif
  28.  
  29. #include <ctime>
  30.  
  31. #include <string>
  32.  
  33. #include "f77-fcn.h"
  34.  
  35. #include "defun-dld.h"
  36. #include "error.h"
  37. #include "gripes.h"
  38. #include "help.h"
  39. #include "mappers.h"
  40. #include "oct-obj.h"
  41. #include "unwind-prot.h"
  42. #include "utils.h"
  43.  
  44. // Possible distributions of random numbers.  This was handled with an
  45. // enum, but unwind_protecting that doesn't work so well.
  46. #define uniform_dist 1
  47. #define normal_dist 2
  48.  
  49. // Current distribution of random numbers.
  50. static int current_distribution = uniform_dist;
  51.  
  52. // Has the seed been set yet?
  53. static int initialized = 0;
  54.  
  55. extern "C"
  56. {
  57.   int *F77_FCN (dgennor, DGENNOR) (const double&, const double&,
  58.                    double&);
  59.  
  60.   int *F77_FCN (dgenunf, DGENUNF) (const double&, const double&,
  61.                    double&);
  62.  
  63.   int *F77_FCN (setall, SETALL) (const int&, const int&);
  64.  
  65.   int *F77_FCN (getsd, GETSD) (int&, int&);
  66.  
  67.   int *F77_FCN (setsd, SETSD) (const int&, const int&);
  68.  
  69.   int *F77_FCN (setcgn, SETCGN) (const int&);
  70. }
  71.  
  72. static double
  73. curr_rand_seed (void)
  74. {
  75.   union d2i { double d; int i[2]; };
  76.   union d2i u;
  77.   F77_XFCN (getsd, GETSD, (u.i[0], u.i[1]));
  78.   return u.d;
  79. }
  80.  
  81. static int
  82. force_to_fit_range (int i, int lo, int hi)
  83. {
  84.   assert (hi > lo && lo >= 0 && hi > lo);
  85.  
  86.   i = i > 0 ? i : -i;
  87.  
  88.   if (i < lo)
  89.     i = lo;
  90.   else if (i > hi)
  91.     i = i % hi;
  92.  
  93.   return i;
  94. }
  95.  
  96. static void
  97. set_rand_seed (double val)
  98. {
  99.   union d2i { double d; int i[2]; };
  100.   union d2i u;
  101.   u.d = val;
  102.   int i0 = force_to_fit_range (u.i[0], 1, 2147483563);
  103.   int i1 = force_to_fit_range (u.i[1], 1, 2147483399);
  104.   F77_XFCN (setsd, SETSD, (i0, i1));
  105. }
  106.  
  107. static char *
  108. curr_rand_dist (void)
  109. {
  110.   if (current_distribution == uniform_dist)
  111.     return "uniform";
  112.   else if (current_distribution == normal_dist)
  113.     return "normal";
  114.   else
  115.     {
  116.       panic_impossible ();
  117.       return 0;
  118.     }
  119. }
  120.  
  121. // Make the random number generator give us a different sequence every
  122. // time we start octave unless we specifically set the seed.  The
  123. // technique used below will cycle monthly, but it it does seem to
  124. // work ok to give fairly different seeds each time Octave starts.
  125.  
  126. static void
  127. do_initialization (void)
  128. {
  129.   time_t now;
  130.   struct tm *tm;
  131.  
  132.   time (&now);
  133.   tm = localtime (&now);
  134.  
  135.   int hour = tm->tm_hour + 1;
  136.   int minute = tm->tm_min + 1;
  137.   int second = tm->tm_sec + 1;
  138.  
  139.   int s0 = tm->tm_mday * hour * minute * second;
  140.   int s1 = hour * minute * second;
  141.  
  142.   s0 = force_to_fit_range (s0, 1, 2147483563);
  143.   s1 = force_to_fit_range (s1, 1, 2147483399);
  144.  
  145.   F77_XFCN (setall, SETALL, (s0, s1));
  146.  
  147.   initialized = 1;
  148. }
  149.  
  150. static octave_value_list
  151. do_rand (const octave_value_list& args, int nargin)
  152. {
  153.   octave_value_list retval;
  154.  
  155.   int n = 0;
  156.   int m = 0;
  157.  
  158.   if (nargin == 0)
  159.     {
  160.       n = 1;
  161.       m = 1;
  162.  
  163.       goto gen_matrix;
  164.     }
  165.   else if (nargin == 1)
  166.     {
  167.       octave_value tmp = args(0);
  168.  
  169.       if (tmp.is_string ())
  170.     {
  171.       string s_arg = tmp.string_value ();
  172.  
  173.       if (s_arg == "dist")
  174.         {
  175.           retval(0) = curr_rand_dist ();
  176.         }
  177.       else if (s_arg == "seed")
  178.         {
  179.           retval(0) = curr_rand_seed ();
  180.         }
  181.       else if (s_arg == "uniform")
  182.         {
  183.           current_distribution = uniform_dist;
  184.  
  185.           F77_XFCN (setcgn, SETCGN, (uniform_dist));
  186.         }
  187.       else if (s_arg == "normal")
  188.         {
  189.           current_distribution = normal_dist;
  190.  
  191.           F77_XFCN (setcgn, SETCGN, (normal_dist));
  192.         }
  193.       else
  194.         error ("rand: unrecognized string argument");
  195.     }
  196.       else if (tmp.is_scalar_type ())
  197.     {
  198.       double dval = tmp.double_value ();
  199.  
  200.       if (xisnan (dval))
  201.         {
  202.           error ("rand: NaN is invalid a matrix dimension");
  203.         }
  204.       else
  205.         {
  206.           m = n = NINT (tmp.double_value ());
  207.  
  208.           if (! error_state)
  209.         goto gen_matrix;
  210.         }
  211.     }
  212.       else if (tmp.is_range ())
  213.     {
  214.       Range r = tmp.range_value ();
  215.       n = 1;
  216.       m = r.nelem ();
  217.       goto gen_matrix;
  218.     }
  219.       else if (tmp.is_matrix_type ())
  220.     {
  221.       // XXX FIXME XXX -- this should probably use the function
  222.       // from data.cc.
  223.  
  224.       Matrix a = args(0).matrix_value ();
  225.  
  226.       if (error_state)
  227.         return retval;
  228.  
  229.       n = a.rows ();
  230.       m = a.columns ();
  231.  
  232.       if (n == 1 && m == 2)
  233.         {
  234.           n = NINT (a (0, 0));
  235.           m = NINT (a (0, 1));
  236.         }
  237.       else if (n == 2 && m == 1)
  238.         {
  239.           n = NINT (a (0, 0));
  240.           m = NINT (a (1, 0));
  241.         }
  242.       else
  243.         warning ("rand (A): use rand (size (A)) instead");
  244.  
  245.       goto gen_matrix;
  246.     }
  247.       else
  248.     {
  249.       gripe_wrong_type_arg ("rand", tmp);
  250.       return retval;
  251.     }
  252.     }
  253.   else if (nargin == 2)
  254.     {
  255.       if (args(0).is_string ())
  256.     {
  257.       if (args(0).string_value () == "seed")
  258.         {
  259.           double d = args(1).double_value ();
  260.  
  261.           if (! error_state)
  262.         set_rand_seed (d);
  263.         }
  264.     }
  265.       else
  266.     {
  267.       double dval = args(0).double_value ();
  268.  
  269.       if (xisnan (dval))
  270.         {
  271.           error ("rand: NaN is invalid as a matrix dimension");
  272.         }
  273.       else
  274.         {
  275.           n = NINT (dval);
  276.  
  277.           if (! error_state)
  278.         {
  279.           m = NINT (args(1).double_value ());
  280.  
  281.           if (! error_state)
  282.             goto gen_matrix;
  283.         }
  284.         }
  285.     }
  286.     }
  287.  
  288.   return retval;
  289.  
  290.  gen_matrix:
  291.  
  292.   if (n == 0 || m == 0)
  293.     {
  294.       Matrix m;
  295.       retval.resize (1, m);
  296.     }
  297.   else if (n > 0 && m > 0)
  298.     {
  299.       Matrix rand_mat (n, m);
  300.       for (int j = 0; j < m; j++)
  301.     for (int i = 0; i < n; i++)
  302.       {
  303.         double val;
  304.         switch (current_distribution)
  305.           {
  306.           case uniform_dist:
  307.         F77_XFCN (dgenunf, DGENUNF, (0.0, 1.0, val));
  308.         rand_mat (i, j) = val;
  309.         break;
  310.  
  311.           case normal_dist:
  312.         F77_XFCN (dgennor, DGENNOR, (0.0, 1.0, val));
  313.         rand_mat (i, j) = val;
  314.         break;
  315.  
  316.           default:
  317.         panic_impossible ();
  318.         break;
  319.           }
  320.       }
  321.  
  322.       retval(0) = rand_mat;
  323.     }
  324.   else
  325.     error ("rand: invalid negative argument");
  326.  
  327.   return retval;
  328. }
  329.  
  330. DEFUN_DLD (rand, args, nargout,
  331.   "rand            -- generate a random value from a uniform distribution\n\
  332. \n\
  333. rand (N)        -- generate N x N matrix\n\
  334. rand (size (A)) -- generate matrix the size of A\n\
  335. rand (N, M)     -- generate N x M matrix\n\
  336. rand (SEED)     -- get current seed\n\
  337. rand (SEED, N)  -- set seed\n\
  338. \n\
  339. See also: randn")
  340. {
  341.   octave_value_list retval;
  342.  
  343.   int nargin = args.length ();
  344.  
  345.   if (nargin > 2 || nargout > 1)
  346.     print_usage ("rand");
  347.   else
  348.     {
  349.       if (! initialized)
  350.     do_initialization ();
  351.  
  352.       retval = do_rand (args, nargin);
  353.     }
  354.  
  355.   return retval;
  356. }
  357.  
  358. static void
  359. reset_rand_generator (void *)
  360. {
  361.   F77_XFCN (setcgn, SETCGN, (current_distribution));
  362. }
  363.  
  364. DEFUN_DLD (randn, args, nargout,
  365.   "randn            -- generate a random value from a normal distribution\n\
  366. \n\
  367. randn (N)        -- generate N x N matrix\n\
  368. randn (size (A)) -- generate matrix the size of A\n\
  369. randn (N, M)     -- generate N x M matrix\n\
  370. randn (SEED)     -- get current seed\n\
  371. randn (SEED, N)  -- set seed\n\
  372. \n\
  373. See also: rand")
  374. {
  375.   octave_value_list retval;
  376.  
  377.   int nargin = args.length ();
  378.  
  379.   if (nargin > 2 || nargout > 1)
  380.     print_usage ("randn");
  381.   else
  382.     {
  383.       if (! initialized)
  384.     do_initialization ();
  385.  
  386.       begin_unwind_frame ("randn");
  387.  
  388.       // This relies on the fact that elements are popped from the
  389.       // unwind stack in the reverse of the order they are pushed
  390.       // (i.e. current_distribution will be reset before calling
  391.       // reset_rand_generator()).
  392.  
  393.       add_unwind_protect (reset_rand_generator, 0);
  394.       unwind_protect_int (current_distribution);
  395.  
  396.       current_distribution = normal_dist;
  397.  
  398.       F77_XFCN (setcgn, SETCGN, (normal_dist));
  399.  
  400.       retval = do_rand (args, nargin);
  401.  
  402.       run_unwind_frame ("randn");
  403.     }
  404.  
  405.   return retval;
  406. }
  407.  
  408. /*
  409. ;;; Local Variables: ***
  410. ;;; mode: C++ ***
  411. ;;; End: ***
  412. */
  413.